Take Home Exercise 01

Published

February 3, 2024

Project Overview

In this take-home exercise, I will be finding the geographical and spatio-temporal distribution of Grab hailing service locations in Singapore.

The distribution will be analysed using the following, which will be displayed on the openstreetmap of Singapore:

  • Traditional Kernel Density Estimation (KDE) layers

  • Network Kernel Density Estimation (NKDE) or Temporal Network Kernel Density Estimation (TNKDE)

Loading the packages

In this section, I will be loading the following packages:

  • maptools

  • sf

  • raster

  • spatstat

  • tmap

  • tidyverse

  • arrow

  • lubridate

  • spNetwork

  • classInt

  • viridis

Show the code
pacman::p_load(maptools, sf, raster, spatstat, tmap, tidyverse, arrow, lubridate, spNetwork, classInt)

Traditional Kernel Density Estimation

Extraction of the Coastal Outline of Mainland Singapore

Now, we can use the data for Singapore’s outline. The data can be found on data.gov.sg and is the Master Plan Subzone Boundary 2014.

Show the code
mpsz_sf <- st_read(dsn = "../../data/data",
                layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `D:\KrisLBT\IS415-GAA\data\data' using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

We have to check the CRS to ensure that it is in the correct coordinates system.

Show the code
st_crs(mpsz_sf)

Since it is in the wrong CRS, we have to correct it to SVY21.

Show the code
mpsz_sf <- st_transform(mpsz_sf, 3414)

Now, we can plot the island.

Show the code
plot(mpsz_sf)

We notice that the planning subzone includes areas where Grab does not serve due to lack of bridges. Hence, any roads in unavailable areas would either not be useful or actually distort the network analysis.

Most of the areas unserved by Grab are the outer islands. For this, we can identify which of the outer islands to exclude from our analysis.

Show the code
outer_islands<-  filter(mpsz_sf, grepl('ISLAND', PLN_AREA_N))
plot(outer_islands["SUBZONE_N"])

However, we realise that not all the islands are unavailable to Grab drivers. Grab drivers are able to pick up and drop off at Sentosa.

Since 2019, Grab drivers have been allowed entry to Jurong Island, provided they have a security pass. However, given that very few drivers are allowed entry and the fact that this region is generally inaccessible to the public, we will be excluding them.

Hence, the islands we need to exclude are:

  • Semakau
  • Southern Group
  • Sudong
  • North-eastern Islands
  • Jurong Island and Bukom
Show the code
serviced_area <- mpsz_sf %>%
  filter(!grepl("SEMAKAU|SOUTHERN GROUP|SUDONG|NORTH-EASTERN ISLANDS|JURONG ISLAND AND BUKOM", SUBZONE_N)) 
plot(serviced_area)

Now, we can blend away the boundaries.

Show the code
sg_sf <- serviced_area %>%
  st_union()

Plot the sg_sf and we now see we have the CoastalOutline of Singapore.

Show the code
plot(sg_sf)

We can save it as the CoastalOutline as a shapefile and read it again

Show the code
# writing the CoastalOutline
st_write(sg_sf, 'data/CoastalOutline.shp')
Show the code
# reading the CoastalOutline
sg_sf <- st_read(dsn='data/',
                 layer = "CoastalOutline")
Reading layer `CoastalOutline' from data source 
  `D:\KrisLBT\IS415-GAA\Take_Home_Exercises\Take_Home_Exercise_1\data' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 21494.3 xmax: 55941.94 ymax: 50256.33
Projected CRS: SVY21 / Singapore TM

Handling of the Grab Posisi data

Reading the Grab Posisi Data

I will now be loading the Grab Posisi data

Show the code
df <- read_parquet(file="../../data/data/GrabPosisi/part-00000.snappy.parquet")
df1 <- read_parquet(file="../../data/data/GrabPosisi/part-00001.snappy.parquet")
df2 <- read_parquet(file="../../data/data/GrabPosisi/part-00002.snappy.parquet")
df3 <- read_parquet(file="../../data/data/GrabPosisi/part-00003.snappy.parquet")
df4 <- read_parquet(file="../../data/data/GrabPosisi/part-00004.snappy.parquet")
df5 <- read_parquet(file="../../data/data/GrabPosisi/part-00005.snappy.parquet")
df6 <- read_parquet(file="../../data/data/GrabPosisi/part-00006.snappy.parquet")
df7 <- read_parquet(file="../../data/data/GrabPosisi/part-00007.snappy.parquet")
df8 <- read_parquet(file="../../data/data/GrabPosisi/part-00008.snappy.parquet")
df9 <- read_parquet(file="../../data/data/GrabPosisi/part-00009.snappy.parquet")

Handling the Grab Posisi data

Time stamp

There was supposed to be a datetime stamp but it wasn’t in that format. As such, the data type needs to be changed.

Show the code
df$pingtimestamp <- as_datetime(df$pingtimestamp)
df1$pingtimestamp <- as_datetime(df1$pingtimestamp)
df2$pingtimestamp <- as_datetime(df2$pingtimestamp)
df3$pingtimestamp <- as_datetime(df3$pingtimestamp)
df4$pingtimestamp <- as_datetime(df4$pingtimestamp)
df5$pingtimestamp <- as_datetime(df5$pingtimestamp)
df6$pingtimestamp <- as_datetime(df6$pingtimestamp)
df7$pingtimestamp <- as_datetime(df7$pingtimestamp)
df8$pingtimestamp <- as_datetime(df8$pingtimestamp)
df9$pingtimestamp <- as_datetime(df9$pingtimestamp)

Merging them into 1 data frame

Merging the origin data frames into 1 origin data frame

For easier handling, we can merge the data frames together.

Show the code
origin_df_all <- df %>%
  rbind(df1) %>%
  rbind(df2) %>%
  rbind(df3) %>%
  rbind(df4) %>%
  rbind(df5) %>%
  rbind(df6) %>%
  rbind(df7) %>%
  rbind(df8) %>%
  rbind(df9)

Extraction of origin and destination

Extraction of origin

Now, I will be extracting the origin.

The way I did it below is by performing the following steps

  • arrange it according to the trj_id (unique id determining where the person wants to at that point)
  • arranging it according to the time stamp (beginning and end)
  • getting the first row - mutating the data to weekday, start hour and the day
Show the code
origin_df <- origin_df_all %>%
  group_by(trj_id) %>%
  arrange(pingtimestamp) %>%
  filter(row_number()==1) %>%
  mutate(weekday = wday(pingtimestamp,
                        label=TRUE,
                        abbr=TRUE),
         start_hr= factor(hour(pingtimestamp)),
         day = factor(mday(pingtimestamp)))
Writing of origin data to rds

Now, I will write the origin data to rds for easier future handling.

Show the code
write_rds(origin_df, "data/rds/origin_dfs.rds")

Reading of RDS data

Now, we can clear the data in our environment and just read the rds data we saved.

Reading for origin data

Show the code
origin_df <- read_rds("data/rds/origin_dfs.rds")

Analysis of Origins

Visualising the Frequency Distribution

Show the code
ggplot(data=origin_df,
       aes(x=weekday)) +
  geom_bar()

As can be seen, the number of trips daily seem to be fairly equally distributed throughout the week.

Visualising the Origins as a Point Map Symbol

We can visualise the geospatial distribution of the origin points.

Show the code
origins_sf <- st_as_sf(origin_df, coords = c("rawlng", "rawlat"),
                       crs=4326) %>%
  st_transform(crs=3414)

Now, we can visualise it using the code chunk below:

Show the code
tm_shape(mpsz_sf)+
  tm_polygons()+
tm_shape(origins_sf) +
  tm_dots()

Convert sf data frames to sp’s Spatial class

We can convert these data to sp’s Spatial* class and then convert into a generic sp format.

Show the code
sg <- as_Spatial(sg_sf)
mpsz <- as_Spatial(mpsz_sf)
origins <- as_Spatial(origins_sf) 

And into SpatialPolygon format

Show the code
sg_sp <- as(sg, "SpatialPolygons")
origins_sp <- as(origins, "SpatialPoints")

We then convert it to ppp format.

Show the code
origins_ppp <- as(origins_sp, "ppp")
origins_ppp
Planar point pattern: 28000 points
window: rectangle = [3628.24, 49845.23] x [25198.14, 49689.64] units

Now, I will plot origins_ppp.

Show the code
plot(origins_ppp)

We must check for duplicates in the data using the code chunk below:

Show the code
any(duplicated(origins_ppp))
[1] FALSE

As there are no duplicates within the data, we do not have to apply any more transformation to the data.

Create owin data

Now, we can convert the Coastal Outline to owin data and plot it out.

Show the code
sg_owin <- as(sg_sp, "owin")
plot(sg_owin)

Combining events

We will now extract Grab origins within the Singapore CoastalOutline

Show the code
origins_SG_ppp = origins_ppp[sg_owin]
Show the code
plot(origins_SG_ppp)

Traditional Kernel Analysis

Now we can begin the traditional kernel analysis

Country level analysis

We can begin by trying to find all of the best kernel and bandwidths for the purposes of our analysis. To do this, we need to find the one with the tightest clusters.

In R, we have 4 possible kernel density bandwidths: diggle, ppl, scott and CvL. All of them will be test below:

Show the code
kde_origins_SG.bw <- density(origins_SG_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 
kde_origins_SG.ppl <- density(origins_SG_ppp, 
                               sigma=bw.ppl, 
                               edge=TRUE,
                               kernel="gaussian")
kde_origins_SG.CvL <- density(origins_SG_ppp, 
                               sigma=bw.CvL, 
                               edge=TRUE,
                               kernel="gaussian")
kde_origins_SG.scott <- density(origins_SG_ppp, 
                               sigma=bw.scott, 
                               edge=TRUE,
                               kernel="gaussian")

We may also plot them out:

Show the code
par(mfrow=c(2,2))
plot(kde_origins_SG.bw)
plot(kde_origins_SG.ppl)
plot(kde_origins_SG.CvL)
plot(kde_origins_SG.scott)

Now, we can check for the tightness of these clusters:

Show the code
# tightness of diggle
bw.diggle(origins_SG_ppp)

# tightness of ppl
bw.ppl(origins_SG_ppp)

# tightness of CvL
bw.CvL(origins_SG_ppp)

# tightness of scott
bw.scott(origins_SG_ppp)

We can observe that the smallest sigma value is from diggle, suggesting that it has the tightest clusters of all the bandwidth. Hence, we will continue using diggle for this analysis.

We can also test for different kernels:

Show the code
par(mfrow=c(2,2))
plot(density(origins_SG_ppp, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"), 
     main="Gaussian")
plot(density(origins_SG_ppp, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="epanechnikov"), 
     main="Epanechnikov")
plot(density(origins_SG_ppp, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="quartic"), 
     main="Quartic")
plot(density(origins_SG_ppp, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="disc"), 
     main="Disc")

As they are identical, there is no need to choose a specific one.

Regions of Interest

We can tighten our analysis to include only areas of interest. From the figures above, we observe higher densities in the planning areas Changi, Jurong East, Woodlands and Marine Parade, we can filter out mpsz to get them.

Show the code
# extraction of Changi Airport
CA = mpsz_sf[mpsz_sf$PLN_AREA_N=="CHANGI",]

# extraction of Jurong East
JE = mpsz_sf[mpsz_sf$PLN_AREA_N=="JURONG EAST",]

# extraction of Woodlands
WL = mpsz_sf[mpsz_sf$PLN_AREA_N=="WOODLANDS",]

# extraction of Marine Parade
MP = mpsz_sf[mpsz_sf$PLN_AREA_N=='MARINE PARADE',]

With this, we now have to perform the same functions we did in Create Owin Data.

Show the code
# turn them into spatial
CA_spatial <- as_Spatial(CA)
JE_spatial <- as_Spatial(JE)
WL_spatial <- as_Spatial(WL)
MP_spatial <- as_Spatial(MP)

# turn them into SpatialPolygons
CA_sp <- as(CA_spatial, "SpatialPolygons")
JE_sp <- as(JE_spatial, "SpatialPolygons")
WL_sp <- as(WL_spatial, "SpatialPolygons")
MP_sp <- as(MP_spatial, "SpatialPolygons")

# convert to owin

CA_owin = as(CA_sp, "owin")
JE_owin = as(JE_sp, "owin")
WL_owin = as(WL_sp, "owin")
MP_owin = as(MP_sp, "owin")

From here, we can extract the events that occurred in these planning areas.

Show the code
origins_CA_ppp = origins_SG_ppp[CA_owin]
origins_JE_ppp = origins_SG_ppp[JE_owin]
origins_WL_ppp = origins_SG_ppp[WL_owin]
origins_MP_ppp = origins_SG_ppp[MP_owin]

We can now perform the same analysis that we did above:

Show the code
kde_origins_CA <- density(origins_CA_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 
kde_origins_JE <- density(origins_JE_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 
kde_origins_WL <- density(origins_WL_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 
kde_origins_MP <- density(origins_MP_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 

We can also plot them out:

Show the code
par(mfrow=c(2,2))
plot(kde_origins_CA,
     main="Changi")
plot(kde_origins_WL,
     main="Woodlands")
plot(kde_origins_JE,
     main="Jurong East")
plot(kde_origins_MP,
     main="Marine Parade")

Clark and Evans test

Entire Singapore

Show the code
clarkevans.test(origins_SG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("clustered"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Z-test

data:  origins_SG_ppp
R = 0.28039, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

Changi

Show the code
clarkevans.test(origins_CA_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Z-test

data:  origins_CA_ppp
R = 0.13547, p-value < 2.2e-16
alternative hypothesis: two-sided

Woodlands

Show the code
clarkevans.test(origins_WL_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Z-test

data:  origins_WL_ppp
R = 0.31778, p-value < 2.2e-16
alternative hypothesis: two-sided

Jurong East

Show the code
clarkevans.test(origins_JE_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Z-test

data:  origins_JE_ppp
R = 0.25797, p-value < 2.2e-16
alternative hypothesis: two-sided

Marine Parade

Show the code
clarkevans.test(origins_MP_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Z-test

data:  origins_MP_ppp
R = 0.51201, p-value < 2.2e-16
alternative hypothesis: two-sided

Network Kernel Density Estimation

While the above does show the distribution of the points and which areas are more concentrated than others, they do not account for the network of the roads. In this section, we will be finding the networks (roads) that are more densely used as origins points in Changi Airport and Marina East.

Data Handling

Handling of Road Data

Now, we can read the road data as provided by openstreetmap. Note that it provides data for Malaysia, Singapore and Brunei all at once.

Show the code
all_roads <- st_read(dsn='../../data/data/data',
                     layer = 'gis_osm_roads_free_1')

We realise that it is in WGS 84, not in SVY21 so we have to transform the data.

Show the code
all_roads <- st_transform(all_roads, 3414)

We can check the CRS again for all_roads.

Show the code
st_crs(all_roads)

Since we’re only interested in the roads in mainland Singapore, we have to filter the roads.

Show the code
SG_roads <- st_intersection(all_roads, sg_sf)
Show the code
st_write(SG_roads, 'data/SG_roads.shp')

We can now read it from here.

Show the code
SG_roads <- st_read(dsn='data', layer= 'SG_roads')
Reading layer `SG_roads' from data source 
  `D:\KrisLBT\IS415-GAA\Take_Home_Exercises\Take_Home_Exercise_1\data' 
  using driver `ESRI Shapefile'
Simple feature collection with 227969 features and 10 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 2679.373 ymin: 23099.51 xmax: 50957.8 ymax: 50220.06
Projected CRS: SVY21 / Singapore TM

Getting the Roads within the Subzones

Now, we can find the streets that exist only inside the subzones.

Show the code
CA_roads <- st_intersection(SG_roads, CA)
WL_roads <- st_intersection(SG_roads, WL)
JE_roads <- st_intersection(SG_roads, JE)
MP_roads <- st_intersection(SG_roads, MP)
Show the code
CA_roads
Simple feature collection with 4007 features and 25 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 42577.65 ymin: 32548.05 xmax: 50286.8 ymax: 41649.2
Projected CRS: SVY21 / Singapore TM
First 10 features:
       osm_id code      fclass                  name  ref oneway maxspeed layer
356  22617051 5113     primary         Loyang Avenue <NA>      F       70     0
360  22617159 5113     primary         Loyang Avenue <NA>      F       70     0
1348 22981700 5113     primary       Telok Paku Road <NA>      F       50     0
3201 34403286 5114   secondary            Loyang Way <NA>      F       50     0
3207 34403387 5122 residential Changi North Street 1 <NA>      B        0     0
3208 34403405 5122 residential Changi North Crescent <NA>      B       50     0
3209 34403417 5122 residential     Changi North Rise <NA>      B        0     0
4001 42063713 5141     service                  <NA> <NA>      F        0     0
4002 42063715 5141     service                  <NA> <NA>      B        0     0
4003 42063716 5141     service                  <NA> <NA>      B        0     0
     bridge tunnel OBJECTID SUBZONE_NO   SUBZONE_N SUBZONE_C CA_IND PLN_AREA_N
356       F      F      221          2 CHANGI WEST    CHSZ02      N     CHANGI
360       F      F      221          2 CHANGI WEST    CHSZ02      N     CHANGI
1348      F      F      221          2 CHANGI WEST    CHSZ02      N     CHANGI
3201      F      F      221          2 CHANGI WEST    CHSZ02      N     CHANGI
3207      F      F      221          2 CHANGI WEST    CHSZ02      N     CHANGI
3208      F      F      221          2 CHANGI WEST    CHSZ02      N     CHANGI
3209      F      F      221          2 CHANGI WEST    CHSZ02      N     CHANGI
4001      F      F      221          2 CHANGI WEST    CHSZ02      N     CHANGI
4002      F      F      221          2 CHANGI WEST    CHSZ02      N     CHANGI
4003      F      F      221          2 CHANGI WEST    CHSZ02      N     CHANGI
     PLN_AREA_C    REGION_N REGION_C          INC_CRC FMEL_UPD_D   X_ADDR
356          CH EAST REGION       ER 7460F2CFB1D7D36C 2014-12-05 44092.34
360          CH EAST REGION       ER 7460F2CFB1D7D36C 2014-12-05 44092.34
1348         CH EAST REGION       ER 7460F2CFB1D7D36C 2014-12-05 44092.34
3201         CH EAST REGION       ER 7460F2CFB1D7D36C 2014-12-05 44092.34
3207         CH EAST REGION       ER 7460F2CFB1D7D36C 2014-12-05 44092.34
3208         CH EAST REGION       ER 7460F2CFB1D7D36C 2014-12-05 44092.34
3209         CH EAST REGION       ER 7460F2CFB1D7D36C 2014-12-05 44092.34
4001         CH EAST REGION       ER 7460F2CFB1D7D36C 2014-12-05 44092.34
4002         CH EAST REGION       ER 7460F2CFB1D7D36C 2014-12-05 44092.34
4003         CH EAST REGION       ER 7460F2CFB1D7D36C 2014-12-05 44092.34
       Y_ADDR SHAPE_Leng SHAPE_Area                       geometry
356  38928.66   14918.11    4848517 LINESTRING (44617.88 41088....
360  38928.66   14918.11    4848517 LINESTRING (43851.2 39632.9...
1348 38928.66   14918.11    4848517 LINESTRING (45352.67 41113....
3201 38928.66   14918.11    4848517 LINESTRING (43761.28 39560....
3207 38928.66   14918.11    4848517 LINESTRING (43152.23 36906....
3208 38928.66   14918.11    4848517 LINESTRING (42843.89 36842....
3209 38928.66   14918.11    4848517 LINESTRING (43130.89 36812....
4001 38928.66   14918.11    4848517 LINESTRING (43433.23 37663....
4002 38928.66   14918.11    4848517 LINESTRING (43436.11 37657....
4003 38928.66   14918.11    4848517 LINESTRING (43511.5 37906.7...
Show the code
WL_roads
Simple feature collection with 7165 features and 25 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 20613.4 ymin: 44814.83 xmax: 25593.19 ymax: 49200.77
Projected CRS: SVY21 / Singapore TM
First 10 features:
       osm_id code        fclass                name  ref oneway maxspeed layer
796  22773625 5115      tertiary  Woodlands Avenue 6 <NA>      F       40     0
819  22774351 5131 motorway_link                <NA> <NA>      F       50     0
825  22774532 5122   residential  Woodlands Drive 14 <NA>      F       40     0
826  22774537 5122   residential  Woodlands Drive 53 <NA>      B       50     0
827  22774541 5122   residential  Woodlands Drive 43 <NA>      F       40     0
832  22775173 5115      tertiary  Woodlands Avenue 5 <NA>      F       60     0
834  22775385 5113       primary Woodlands Avenue 12 <NA>      F       70     0
835  22775386 5114     secondary  Woodlands Avenue 1 <NA>      F       50     0
836  22775389 5114     secondary  Woodlands Avenue 1 <NA>      F       50     0
3703 37584630 5111      motorway  Seletar Expressway  SLE      F       90     1
     bridge tunnel OBJECTID SUBZONE_NO       SUBZONE_N SUBZONE_C CA_IND
796       F      F      282          4 WOODLANDS SOUTH    WDSZ04      N
819       F      F      282          4 WOODLANDS SOUTH    WDSZ04      N
825       F      F      282          4 WOODLANDS SOUTH    WDSZ04      N
826       F      F      282          4 WOODLANDS SOUTH    WDSZ04      N
827       F      F      282          4 WOODLANDS SOUTH    WDSZ04      N
832       F      F      282          4 WOODLANDS SOUTH    WDSZ04      N
834       F      F      282          4 WOODLANDS SOUTH    WDSZ04      N
835       F      F      282          4 WOODLANDS SOUTH    WDSZ04      N
836       F      F      282          4 WOODLANDS SOUTH    WDSZ04      N
3703      T      F      282          4 WOODLANDS SOUTH    WDSZ04      N
     PLN_AREA_N PLN_AREA_C     REGION_N REGION_C          INC_CRC FMEL_UPD_D
796   WOODLANDS         WD NORTH REGION       NR 8A4E14DAC4ACE11C 2014-12-05
819   WOODLANDS         WD NORTH REGION       NR 8A4E14DAC4ACE11C 2014-12-05
825   WOODLANDS         WD NORTH REGION       NR 8A4E14DAC4ACE11C 2014-12-05
826   WOODLANDS         WD NORTH REGION       NR 8A4E14DAC4ACE11C 2014-12-05
827   WOODLANDS         WD NORTH REGION       NR 8A4E14DAC4ACE11C 2014-12-05
832   WOODLANDS         WD NORTH REGION       NR 8A4E14DAC4ACE11C 2014-12-05
834   WOODLANDS         WD NORTH REGION       NR 8A4E14DAC4ACE11C 2014-12-05
835   WOODLANDS         WD NORTH REGION       NR 8A4E14DAC4ACE11C 2014-12-05
836   WOODLANDS         WD NORTH REGION       NR 8A4E14DAC4ACE11C 2014-12-05
3703  WOODLANDS         WD NORTH REGION       NR 8A4E14DAC4ACE11C 2014-12-05
       X_ADDR   Y_ADDR SHAPE_Leng SHAPE_Area                       geometry
796  23609.57 45692.88   5211.384    1576001 LINESTRING (24007.3 45538.1...
819  23609.57 45692.88   5211.384    1576001 LINESTRING (22821.47 45515....
825  23609.57 45692.88   5211.384    1576001 LINESTRING (23444.28 46022....
826  23609.57 45692.88   5211.384    1576001 LINESTRING (23990.45 46088....
827  23609.57 45692.88   5211.384    1576001 LINESTRING (23447.47 46014....
832  23609.57 45692.88   5211.384    1576001 LINESTRING (24645.58 46021....
834  23609.57 45692.88   5211.384    1576001 LINESTRING (23836.85 45038....
835  23609.57 45692.88   5211.384    1576001 LINESTRING (23124.88 45989....
836  23609.57 45692.88   5211.384    1576001 LINESTRING (24156.4 45407.0...
3703 23609.57 45692.88   5211.384    1576001 LINESTRING (23684.88 44959....
Show the code
JE_roads
Simple feature collection with 7131 features and 25 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 14254.68 ymin: 30994.22 xmax: 19398.25 ymax: 37289.81
Projected CRS: SVY21 / Singapore TM
First 10 features:
       osm_id code       fclass                     name  ref oneway maxspeed
1084 22903885 5114    secondary             Penjuru Road <NA>      F       60
1085 22903886 5121 unclassified             Penjuru Road <NA>      F       60
1443 23101520 5112        trunk              Jalan Buroh <NA>      F       70
3791 39959160 5112        trunk              Jalan Buroh <NA>      F       70
3792 39959161 5112        trunk              Jalan Buroh <NA>      F       70
5003 70583774 5115     tertiary              Pandan Road <NA>      B       50
5004 70583780 5121 unclassified             Penjuru Lane <NA>      B       50
5005 70583783 5114    secondary          Tanjong Penjuru <NA>      B       50
5006 70583786 5114    secondary          Tanjong Penjuru <NA>      F       50
5007 70583790 5121 unclassified Tanjong Penjuru Crescent <NA>      B       50
     layer bridge tunnel OBJECTID SUBZONE_NO        SUBZONE_N SUBZONE_C CA_IND
1084     0      F      F       80         10 PENJURU CRESCENT    JESZ10      N
1085     0      F      F       80         10 PENJURU CRESCENT    JESZ10      N
1443     1      T      F       80         10 PENJURU CRESCENT    JESZ10      N
3791     1      T      F       80         10 PENJURU CRESCENT    JESZ10      N
3792     0      F      F       80         10 PENJURU CRESCENT    JESZ10      N
5003     0      F      F       80         10 PENJURU CRESCENT    JESZ10      N
5004     0      F      F       80         10 PENJURU CRESCENT    JESZ10      N
5005     0      F      F       80         10 PENJURU CRESCENT    JESZ10      N
5006     0      F      F       80         10 PENJURU CRESCENT    JESZ10      N
5007     0      F      F       80         10 PENJURU CRESCENT    JESZ10      N
      PLN_AREA_N PLN_AREA_C    REGION_N REGION_C          INC_CRC FMEL_UPD_D
1084 JURONG EAST         JE WEST REGION       WR BC263278706376DE 2014-12-05
1085 JURONG EAST         JE WEST REGION       WR BC263278706376DE 2014-12-05
1443 JURONG EAST         JE WEST REGION       WR BC263278706376DE 2014-12-05
3791 JURONG EAST         JE WEST REGION       WR BC263278706376DE 2014-12-05
3792 JURONG EAST         JE WEST REGION       WR BC263278706376DE 2014-12-05
5003 JURONG EAST         JE WEST REGION       WR BC263278706376DE 2014-12-05
5004 JURONG EAST         JE WEST REGION       WR BC263278706376DE 2014-12-05
5005 JURONG EAST         JE WEST REGION       WR BC263278706376DE 2014-12-05
5006 JURONG EAST         JE WEST REGION       WR BC263278706376DE 2014-12-05
5007 JURONG EAST         JE WEST REGION       WR BC263278706376DE 2014-12-05
       X_ADDR   Y_ADDR SHAPE_Leng SHAPE_Area                       geometry
1084 17651.31 31799.61   9876.049    3049720 LINESTRING (16991.44 32621....
1085 17651.31 31799.61   9876.049    3049720 LINESTRING (16738.48 31589....
1443 17651.31 31799.61   9876.049    3049720 LINESTRING (16290.72 32753....
3791 17651.31 31799.61   9876.049    3049720 LINESTRING (19056.98 32186....
3792 17651.31 31799.61   9876.049    3049720 LINESTRING (18978.73 32113....
5003 17651.31 31799.61   9876.049    3049720 LINESTRING (18411.82 31559....
5004 17651.31 31799.61   9876.049    3049720 LINESTRING (17208.01 32320....
5005 17651.31 31799.61   9876.049    3049720 LINESTRING (16869.64 32131....
5006 17651.31 31799.61   9876.049    3049720 LINESTRING (17506.95 32049....
5007 17651.31 31799.61   9876.049    3049720 LINESTRING (17697.91 31764....
Show the code
MP_roads
Simple feature collection with 2864 features and 25 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 32697.12 ymin: 29763 xmax: 37589.95 ymax: 32843.37
Projected CRS: SVY21 / Singapore TM
First 10 features:
       osm_id code       fclass                         name  ref oneway
922  22810641 5114    secondary    Tanjong Katong Road South <NA>      F
1202 22930502 5113      primary            Marina East Drive <NA>      F
3760 39477475 5111     motorway           East Coast Parkway  ECP      F
3761 39477476 5111     motorway           East Coast Parkway  ECP      F
4078 44122367 5152     cycleway               Park Connector <NA>      B
4097 44133236 5153      footway      Underpass to Meyer Road <NA>      B
4140 44488255 5153      footway                  Katong Park <NA>      B
4141 44488257 5121 unclassified East Coast Park Service Road <NA>      B
4142 44488265 5114    secondary    Tanjong Katong Road South <NA>      F
5509 74729034 5111     motorway           East Coast Parkway  ECP      F
     maxspeed layer bridge tunnel OBJECTID SUBZONE_NO        SUBZONE_N
922        50     0      F      F       44          5 MARINA EAST (MP)
1202       60     0      F      F       44          5 MARINA EAST (MP)
3760       90     1      T      F       44          5 MARINA EAST (MP)
3761       80     0      F      F       44          5 MARINA EAST (MP)
4078        0     0      F      F       44          5 MARINA EAST (MP)
4097        0    -1      F      T       44          5 MARINA EAST (MP)
4140        0    -1      F      T       44          5 MARINA EAST (MP)
4141       50     0      F      F       44          5 MARINA EAST (MP)
4142       50     1      T      F       44          5 MARINA EAST (MP)
5509       90     0      F      F       44          5 MARINA EAST (MP)
     SUBZONE_C CA_IND    PLN_AREA_N PLN_AREA_C       REGION_N REGION_C
922     MPSZ05      N MARINE PARADE         MP CENTRAL REGION       CR
1202    MPSZ05      N MARINE PARADE         MP CENTRAL REGION       CR
3760    MPSZ05      N MARINE PARADE         MP CENTRAL REGION       CR
3761    MPSZ05      N MARINE PARADE         MP CENTRAL REGION       CR
4078    MPSZ05      N MARINE PARADE         MP CENTRAL REGION       CR
4097    MPSZ05      N MARINE PARADE         MP CENTRAL REGION       CR
4140    MPSZ05      N MARINE PARADE         MP CENTRAL REGION       CR
4141    MPSZ05      N MARINE PARADE         MP CENTRAL REGION       CR
4142    MPSZ05      N MARINE PARADE         MP CENTRAL REGION       CR
5509    MPSZ05      N MARINE PARADE         MP CENTRAL REGION       CR
              INC_CRC FMEL_UPD_D  X_ADDR   Y_ADDR SHAPE_Leng SHAPE_Area
922  1575868DF0D30E32 2014-12-05 33715.7 30512.25   6657.151    1590340
1202 1575868DF0D30E32 2014-12-05 33715.7 30512.25   6657.151    1590340
3760 1575868DF0D30E32 2014-12-05 33715.7 30512.25   6657.151    1590340
3761 1575868DF0D30E32 2014-12-05 33715.7 30512.25   6657.151    1590340
4078 1575868DF0D30E32 2014-12-05 33715.7 30512.25   6657.151    1590340
4097 1575868DF0D30E32 2014-12-05 33715.7 30512.25   6657.151    1590340
4140 1575868DF0D30E32 2014-12-05 33715.7 30512.25   6657.151    1590340
4141 1575868DF0D30E32 2014-12-05 33715.7 30512.25   6657.151    1590340
4142 1575868DF0D30E32 2014-12-05 33715.7 30512.25   6657.151    1590340
5509 1575868DF0D30E32 2014-12-05 33715.7 30512.25   6657.151    1590340
                           geometry
922  LINESTRING (35230.6 31057.3...
1202 LINESTRING (33760.43 30638....
3760 LINESTRING (33860.36 30899....
3761 LINESTRING (33706.71 30894....
4078 LINESTRING (34789.2 30864.6...
4097 LINESTRING (34772.07 30914....
4140 LINESTRING (34005.49 30906....
4141 LINESTRING (35255.73 30985,...
4142 LINESTRING (35227.24 31091....
5509 LINESTRING (34854.12 30932....

We notice that all of the above are geometry type geometry and not a linestring. We can convert them with the following code chunk:

Show the code
CA_roads <- CA_roads %>%
  st_cast("LINESTRING")

WL_roads <- WL_roads %>%
  st_cast("LINESTRING")

JE_roads <- JE_roads %>%
  st_cast("LINESTRING")

MP_roads <- MP_roads %>%
  st_cast("LINESTRING")

Extraction of the Events Within the Subzones

Before we can conduct analysis on the network, we will also need to constrict the events to exclusively the ones that occurred within these two subzones.

Show the code
#extraction of origin events that occurred within Changi and Marine Parade
CA_origins <- st_intersection(origins_sf, CA)
WL_origins <- st_intersection(origins_sf, WL)
JE_origins <- st_intersection(origins_sf, JE)
MP_origins <- st_intersection(origins_sf, MP)

We can now plot this data. In this instance, we can use Changi as an example:

Show the code
tmap_mode('view')
tm_shape(CA_origins)+
  tm_dots() +
  tm_shape(CA_roads) +
  tm_lines()
Show the code
tmap_mode('plot')

Network Constrained KDE (NetKDE) Analysis

Preparing the lixel objects

We can now prepare the lixel objects.

Show the code
CA_lixels <- lixelize_lines(CA_roads, 
                         750, 
                         mindist = 375)
WL_lixels <- lixelize_lines(WL_roads, 
                         750, 
                         mindist = 375)
JE_lixels <- lixelize_lines(JE_roads, 
                         750, 
                         mindist = 375)
MP_lixels <- lixelize_lines(MP_roads, 
                         750, 
                         mindist = 375)

Generating Line Points

Show the code
CA_samples <- lines_center(CA_lixels)
WL_samples <- lines_center(WL_lixels)
JE_samples <- lines_center(JE_lixels)
MP_samples <- lines_center(MP_lixels)

Performing NetKDE

Show the code
CA_densities <- nkde(CA_roads, 
                  events = CA_origins,
                  w = rep(1,nrow(CA_origins)),
                  samples = CA_samples,
                  kernel_name = "quartic",
                  bw = 3.407207, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 10, #we aggregate events within a 10m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)
# care about bw (bandwith) and kernel_name